//============================================================================
// Name        : rscj.cpp
// Author      : Bob
// Version     :0
// Copyright   : Your copyright notice
// Description : your mom
//============================================================================


#include "stdafx.h"
#include <vector>
#include <ctime>


using namespace std;

const long double pi=3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211;
const long double tmin=0;
const long double tmax=500;
const long double imax=3;
const long double deltai=0.01;
const long double deltat=0.001;
const long double ynull=1;
const long double betac=0.1;

long double phitemp=0;
long double utemp=0;

void write(vector<long double> data1,vector<long double> data2,vector<long double> data3, int n);
void write(vector<long double> data1,vector<long double> data2, int n);
long double fkt(long double dc1,long double y,long double yt);
long double U_t(long double dc1, long double Unull);


int main() {

	clock_t begin = clock();


	int n=2*(imax/deltai);

	vector<long double> I;
	vector<long double> U;
	I.resize(n);
	U.resize(n);
	I[0]=0;
	U[0]=U_t(I[0],0.0);
	for(int i=1;i<n/2;i++){
		I[i]=I[i-1]+deltai;
		U[i]=U_t(I[i],U[i-1]);
	}
	for(int i=n/2;i<n;i++){
			I[i]=I[i-1]-deltai;
			U[i]=U_t(I[i],U[i-1]);
		}

	write(I,U,n);

	clock_t end = clock();
	double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
	cout<<elapsed_secs<<endl;
	char i;
	cin>>i;

	return elapsed_secs;
}

long double U_t(long double dc1, long double Unull){
	int n=(tmax-tmin)/deltat;
		vector<long double> t;
		vector<long double> v;
		vector<long double> phi;
		t.resize(n);
		v.resize(n);
		phi.resize(n);
		for(int i=0;i<n;i++){
			t[i]=tmin+deltat*i;


		}

		long double v0,phi0,va,vb,vc,dv0,dva,dvb,dvc,phia,phib,phic,dphi0,dphia,dphib,dphic;

			//v[0]=utemp;
			//phi[0]=phitemp;
			//v[0]=0;
			v[0]=Unull;
			phi[0]=0;
			for(int i=0;i<n-1;i++){
				v0=v[i];
				phi0=phi[i];
				dphi0=v0;
				dv0=fkt(dc1,v0,phi0);
				phia=phi0+deltat/2.0*dphi0;
				va=v0+deltat/2.0*dv0;
				dphia=va;
				dva=fkt(dc1,va,phia);
				phib=phi0+deltat/2.0*dphia;
				vb=v0+deltat/2.0*dva;
				dphib=vb;
				dvb=fkt(dc1,vb,phib);
				phic=phi0+deltat*dphib;
				vc=v0+deltat*dvb;
				dphic=vc;
				dvc=fkt(dc1,vc,phic);
				v[i+1]=v0+deltat/6.0*(dv0+2*(dva+dvb)+dvc);
				phi[i+1]=phi0+deltat/6.0*(dphi0+2*(dphia+dphib)+dphic);
			}

		long double U=0;
		for(int i=0;i<n;i++){
			U=U+v[i];
		}
			U=U/n;
		//write(t,v,phi,n);
			utemp=v[n-1];
			if(phi[n-1]>2.0*pi){
			phitemp=phi[n-1];
			int l=phi[n-1]/(2.0*pi);
			phitemp=phi[n-1]-l*2.0*pi;
			}
			
	return U;

}


void write(vector<long double> data1,vector<long double> data2,vector<long double> data3, int n){
	fstream datei;
	datei.open("data.txt", ios::out);
	for(int i=0;i<n;i++){
		datei<<data1[i]<<" "<<data2[i]<<" "<<data3[i]<<endl;
	}
	datei.close();
}

void write(vector<long double> data1,vector<long double> data2, int n){
	fstream datei;
	datei.open("data.txt", ios::out);
	for(int i=0;i<n;i++){
		datei<<data1[i]<<" "<<data2[i]<<endl;
	}
	datei.close();
}



long double fkt(long double dc1,long double y,long double yt){
	return ((1/betac)*(dc1-y-sin(yt)));
}
